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i ABSTRACT 

CN \ We used one- dimensional radiative-transport radiation hydrodynamical sim- 

\ ulations to investigate the formation of clumping in freely-expanding supernova 

ejecta due to the radioactive heating from the 56 Ni — ► 56 Co — > 56 Fe decay 
Tj- | sequence. The heating gives rise to an inflated Nickel bubble, which induces a 

forward shock that compresses the outer ambient gas into a shell. The radiative 
energy deposited by the radioactivity leaks out across the shock by radiative 
diffusion, and we investigate its effect on the evolution of the ejecta structure. 
Compared to the hydrodynamical adiabatic approximation with 7 = 4/3, the 
preshock gas becomes accelerated by the radiation outflow. The shock is thus 
weakened and the shell becomes broader and less dense. The thickness of the 
shell takes up < 4% of the radius of the bubble, and the structure of the shell 
can be approximately described by a self-similar solution. We compared the 
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properties of the shell components with those of the ejecta clumps indicated by 
our previous hydrodynamical simulations for the interaction of clumps with the 
outer supernova remnant. The high density contrast across the shell, x ~ 100? is 
compatible with that of ejecta clumps as indicated for Tycho's knots, but there 
is insufficient dense gas to cause a pronounced protrusion on the outline of a core 
collapse supernova remnant, like the bullets in the Vela remnant. 

Subject headings: radiative transfer — supernova remnants — supernovae: gen- 
eral 



1. INTRODUCTION 



Observations of SN 1987A showed that the distribution of Fe in the ejecta is not what 
would be expected in the simplest models: it extended to higher velocity than expected and 
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had a large filling factor for its mass of 0.07 M determined from the supernova light curve 
(McCray 1993; Li, McCray, & Sunyaev 1993). A plausible mechanism for the large filling 
factor is the Ni bubble effect, in which the radioactive progenitors of the Fe expand relative 
to their surroundings because of the radioactive power deposition (Woosley 1988; Li et al. 
1993; Basko 1994). This effect is important during the first ~ 10 days after the supernova, 
when the radioactive power is significant and the diffusion of energy is not yet an important 
process. 

An expected effect of the Ni bubble expansion is to create clumps in the nonradioactive 
gas. There is widespread evidence that the ejecta of core collapse supernovae are clumpy. The 
oxygen line profiles in the nearby Type II supernovae SN 1987A (Stathakis et al. 1991) and 
SN 1993J (Spyromilio 1994; Matheson et al. 2000) showed evidence for structure, implying 
that the gas is clumped. The velocity range for the emission extends to 1,500 km s -1 in 
SN 1987A and 4, 000 km s" 1 in SN 1993J. Similar evidence for clumping has been found in 
the Type lb supernova SN 1985F (Filippenko & Sargent 1989). Among young supernova 
remnants, Cas A is the prototype of the oxygen-rich SNRs, which show evidence for freely 
expanding, oxygen-rich ejecta in clumps. Puppis A, with an age ~ 3, 700 yr, is a more elderly 
example of such a remnant (Winkler et al. 1988). 

Wang & Chevalier (2001, 2002, hereafter WC01, WC02) investigated the role of clumps 
in the evolution of supernova remnants. In Tycho's remnant, believed to be the remnant of a 
Type la supernova, the presence of ejecta knots near the outer shock front requires a density 
contrast x ^ 100 relative to the surrounding ejecta (WC01). The remarkable protrusions in 
the X-ray image of the Vela remnant are likely to be caused by ejecta clumps (Aschenbach 
et al. 1995) and WC02 found that x ~ 1000 is needed to create the structures. The free 
expansion velocities for the clumps were estimated at ~ 3, 000 km s _1 . 

Our aim here is to investigate whether the Ni bubble effect can cause the clump struc- 
ture apparently needed in the Tycho and Vela remnants. We compute one-dimensional 
hydrodynamical calculations, building on the work of Basko (1994), suitable for core col- 
lapse supernovae. We extend Basko's work by including the effects of radiative diffusion and 
by carefully examining the density structure of the shell swept up by the Ni bubble. We fur- 
ther consider Type la supernovae using the same methodology. We show our computational 
setup and methods in § 2. The solutions for radiation- hydrodynamical effects are given in 
§ 3. We also draw on the analogy to power input in a pulsar bubble to provide insight into 
the shell structure and evolution. In § 4 and § 5, we discuss the radiative effects on the 
Rayleigh- Taylor instabilities of the Ni bubble shell and the inferred ejecta-clump properties. 
Our conclusions are in § 6. 
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2. INITIAL CONDITIONS AND METHODS 

We consider all of the 56 Ni synthesized in a supernova explosion initially resides in an 
isolated sphere in the center of supernova ejecta, with a mass M^i and a density contrast 
lv relative to the surrounding diffuse ejecta. The total radioactive energy released in the 
56 Ni — > 56 Co — >• 56 Fe decay sequence is deposited at a rate q(t) (Basko 1994) 

g (t) - ® Ni CXp ( - ^ | Qco^Vi-t/tco) ~ e*P{-t/tNi) ^ 

tm \ tm J tco 1 — tm/tco 

where tm^Co an d QniiQco are, respectively, the mean life and specific decay energy of 
56 Ni, 56 Co: 

t m = 7.6 x 10 5 sec = 8.8 days f 56 Ni -> 56 Co 

Q Ni = 2.14 MeV/decay = 3.69 x 10 16 ergs g" 1 

t Co = 9.64 x 10 6 sec = 111.5 days 56 Co -»• 56 Fe 

Q Co = 4.57 Mev/decay = 7.87 x 10 16 ergs g" 1 . 

(Lide 1992; f 8.5 days, Fireston & Shirley 1996). The accumulation of the total radioactive 
energy rises from 10 4 sec and saturates at > 10 7 sec (about 25% of the total radiative energy 
is deposited between 10 7 and 10 8 sec, Fig. 1). Before the saturation point is reached, the 
radioactive energy can leak out of the bubble by radiative diffusion. This can occur while 
the gas remains optically thick to the ~ 1 MeV 7-rays. 

We assume that the diffuse supernova ejecta are initially freely-expanding so that each 
gas element moves with a constant velocity v = r/t. The part of the ejecta unaffected by the 
radioactivity has its density drop as t~ 3 , and the pressure drops as t -4 for a 7 = 4/3, radiation 
dominated fluid. We first refer our simulation initial conditions to Basko's parameters for 
SN 1987A (Basko 1994 and the references therein): 

M Ni = 0.075 M Q , tv = 3 

t = 10 4 sec, p a0 = 10~ 4 g cm" 3 , P = 2.5 x 10 11 dyn cm" 2 , 

where to, Pao and P are the initial age, ejecta substrate density (surrounding the bubble), 
and background pressure, respectively. The initial velocity at the bubble edge r = i? is 
given by the initial 56 Ni density contrast and age: 

TT R fSM m \ 1/3 1 7 _ 1/3 , 

U — — = — = 7 x 10 iv 16 cm s . 

t \4irpmJ t 

For iv = 3, R = 5.0 x 10 11 cm. The initial pressure is assumed uniformly distributed inside 
and outside the bubble, independent of the bubble size. The pressure in the Ni bubble is 
determined by the radioactive power input at later times. 
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We note that the above parameters give an ejecta structure comparable to the inner 
flat-density component of the power-law ejecta density model that has a power index n = 8 
(p oc r~ n in the outer parts), an ejecta mass M — 10 M , and an explosion energy E = 10 51 
ergs (Chevalier & Liang 1989). The transition between the flat and power law parts of 
the density profile occurs at v — 3162 km s" 1 for these parameters. This set of explosion 
conditions is suitable for core collapse supernovae and had been used in our two-dimensional 
hydrodynamical simulations for the Vela supernova remnant (WC02). The energy density is 
initially dominated by the energy deposited by the supernova shock. Shigeyama & Nomoto's 
(1990, their Fig. 2) model for SN 1987A indicates an adiabat re = 6.3 x 10 15 (cgs units) 
where P = np 1 ] for t = 10 4 sec, p = 9.4 x 10~ 5 g cm -3 in our density model, and so 
Po = 2.7 x 10 10 dynes cm" 2 . We take Pq = 2.7 x 10 10 dynes cm -2 as the appropriate 
value at 10 4 sec. For to — 100 sec, Pq = 2.7 x 10 18 dynes cm" 2 . The time to set up the 
free expansion phase depends on the radius of the star because of the effects of reverse 
shock waves. For a blue supergiant progenitor star, like that of SN 1987A, the timescale 
is ~ 10 4 sec; for a red supergiant progenitor, the timescale is ~ 10 5 sec. It is needed to 
initiate the simulations at t ~ 10 2 sec to acquire a numerically saturated solution at t > 10 4 
sec (the age of interest). The use of a freely-expanding velocity prior to this stage should 
not affect the 1-D symmetrical solution, considering that the radioactive heating is strong. 
Models of supernova explosions predict smooth 56Ni profiles in both thermonuclear and core 
collapse supernovae, with the density contrast u below a factor of 3 (Nomoto, Thielemann 
& Yokoi 1984; Woosley 1988). We vary uj in the range of 0.1 < uo < 3, with the other initial 
parameters scaled accordingly to the standard values described above. 

Next we consider the initial parameters for Type la supernovae. We assume a powerlaw 
density profile n = 8 as in core collapse supernovae, with M — 1.4 Mq, E — 10 51 ergs, and 
Mpji = 0.5 M Q . The transition between the flat and power law parts of the density profile 
occurs at v — 8452 km s" 1 for these parameters. The background pressure in Type la ejecta 
is not well documented. To estimate, we make use of the entropy change dS = dQ/T during 
the nuclear burning from C/O into Ni. During the process, kT ~ 0.6 MeV, and the energy 
released per nucleon is 

dQ ~ [2m(C) + 2m(0) - m(Ni)] c 2 /56 ~ 0.79 MeV/nucleon, (2) 

where m(X) is the mass of element X; m(C) = 12.0000 amu, m(O) = 15.9949 amu, and 
m (Ni) = 55.9421 amu (Audi et al. 2003), and c is the speed of light. Then dS - 
1.3A; /nucleon = 1.8 x 10" 16 erg/K/nucleon, where k is the Boltzmann constant (Shigeyama, 
private communication). The adiabat re was obtained by relating the entropy per nucleon 
S = [(4aT 3 /3)(m(Ni)/56)/p], where the initial entropy before the burning is ignored, to the 
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radiative pressure P = aT 4 /3 and P — npt. Subsituting P and k into S, we find 



K 



'3\l/4/ ^ • 

V l 4(m(Ni)/56)- 



4/3 

~ 6.1 x 10 14 (cgs units). (3) 



For to = 10 4 sec, po — 6.92 x 10~ 7 g cm -3 in our density model, and so P = 3.7 x 10 6 dynes 
cm" 2 . The adiabat (and so the background pressure) could be higher, considering there is 
an initial entropy distribution that acts to freeze the deflagration front. The background 
energy density has the profile of the mass density because of the assumed uniformity. 

We first approach the problem using hydrodynamical simulations with 7 = 4/3 in 
both the bubble interior and the ejecta substrate. We follow the bubble-shell interface on a 
spherical expanding grid, and add the specific radioactive energy to the internal energy within 
the bubble uniformly and locally, i.e., in proportion to the local mass density. A reflecting 
condition is used for the inner boundary, and a non-zero gradient outflow condition is invoked 
on the outer boundary to eliminate the spurious shock raised by the grid expansion. 

We further include radiation hydrodynamics (RHD) in the simulations using the ZEUS2D 
code (Stone, Mihalas, & Norman 1992). The radioactive power is deposited into the radiation 
instead of the material and 7 = 5/3 for the matter. The code is based on finite difference and 
finite volume on an Eulerian grid and uses an artificial viscosity to smooth shock transitions. 
The RHD algorithm employs nonrelativistic full radiative transport, under the assumption 
of LTE (local thermal equilibrium) and gray opacity, and is thus equally applicable to both 
optically thin and thick media. The flow evolves towards the optically thin regime where 
radiation tends to stream away across the Ni bubble shell. The full transport scheme ensures 
the proper evolution while the radioactive energy is continually added within the bubble. 

At the age of interest, the opacity is dominated by electron scattering. A maximum 
absorption 1/1000 times the magnitude of scattering is applied in our calculation. For com- 
pletely ionized heavy elements, the Thompson scattering gives an opacity 0.2 p(r,t) cm -1 , 
where p(r, t) is the material density. The initial energy is equally partitioned between radi- 
ation and material. Considering that a Type la has lower temperature and smaller radius, 
Prad ~ r ~ A an d P gas ~ r~ 5 , the equi-partition assumption may be more valid for core collapse 
super novae. 

In the RHD algorithm, the coupled radiation and material energy densities are solved 
using Newton-Raphson iteration for convergence. When the radiation pressure is in excess 
of the gas pressure by several orders of magnitude, convergence becomes extremely difficult 
to reach. We thus stop the simulations prior to the optically thin stage. In our situation, the 
radiation effects are most important when the gas is still optically thick, so that the diffusion 
approximation is adequate. The choice of the flux limiters from which the Eddington tensor 
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is computed does not affect the diffusion result. For details of radiation hydrodynamics, we 
refer readers to Stone et al. (1992) and the references therein. 



3. SHELL STRUCTURE AND EVOLUTION 
3.1. Simulations 

We show in Figs. 2, 3, and 4 the ejecta structure and evolution properties in the core 
collapse supernova model using n = 8, E = 10 51 erg s _1 , and M = 10 M & . These simulations 
were initiated with t — 10 2 sec, M^i = 0.075 M and ou — 3. The diffuse ejecta have an 
initial density po = 94.4 g cm -3 and an accumulated radioactive energy density e ra do = 
1.38 x 10 15 erg cm -3 . The initial material and radiation energy density in the background 
are equally divided, e = e r0 =1.35 x 10 10 erg cm -3 . 

We first examined the evolution of the flow in the purely hydrodynamic simulations. 
The expansion of the bubble gives rise to a strong forward shock behind which the ambient 
gas is compressed into a dense shell (Figs. 2 and 3). The inner edge of the shell is a contact 
discontinuity where the gas comoves with the bubble-shell interface; it has a high density 
because of the acceleration of the shell. The shock front is continually accelerated outward 
(Fig. 4(a)). To describe the acceleration, we use a = (dR s h/dt)/(R s h/t), the expansion rate 
evaluated at the shock front R s h, equivalent to the velocity ratio of the shell to the preshock 
diffuse ejecta. We find that a rises to a maximum < 1.20 around < 10 6 sec, and subsequently 
tends to 1.0 at late times as the shell becomes comoving with the preshock gas (Fig. 4(b)). 
Because the background pressure has little influence on the dynamics, the shell structure 
becomes frozen into the flow. 

The shell is very thin, with a thickness as a fraction of the radius (3 = h s h/R s h ^$ 0.02, 
where R s h and h s h are respectively the radius and thickness of the shell (Fig. 4(c)). The 
initial shock compression is set by the jump condition 

Pi (7 + l)-M 2 U) 
po ( 7 -l)X 2 + 2' {) 

where 

~ c ~ V dt t ) \ 1Po ) 1 ' 

is the Mach number of the shock, p 1 is the postshock density, v s h is the shock velocity 
relative to the preshock gas, and Co is the sound speed in the preshock medium. For a 
radiation dominated, strong shock, pi/po = 7. The accelerating shell has a higher density in 
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the postshock region so that a typical compression factor in the shell (relative to the preshock 
gas) is x ~ 10 (Fig- 4(d)). The compression factor in the bubble is x ~ 0.1 (Fig. 4(e)). We 
note that (3 ~ (3x) _1 - The highest density occurs at the inner edge, where the density is 
limited by the numerical resolution (Fig. 3). 

In the RHD case, the energy from the heated gas leaks out from the bubble energy 
reservoir as time progresses (Fig. 2). The gas is radiation dominated. Radiative diffusion 
allows the internal energy to propagate ahead of the shock wave, and so eliminates the 
temperature and the radiation pressure gradients. In order for diffusion to be dynamically 
important, it requires that the diffusion time scale, tj, be smaller than the hydrodynamic 
time scale of the flow, th] i.e., 

**f « ^ < « — (6) 
Ac v sh 

where A is the photon mean free path. We note that 



t 



Mt)(t)- < 7 > 



At t « 3 x 10 6 s, h sh = 10 13 cm, A = l/np = 5 x 10 10 cm, so that t d = h 2 sh /Xc = 6 x 10 4 
s. At this time, v s h ~ 1 x 10 7 cm s _1 , so ta/th ~ 0.07. The gas is optically thick, but 
radiative diffusion is important. As a result of diffusion, the preshock gas is accelerated by 
the radiation; the shock is then weakened. The shock front develops a radiative precurso 
(Fig. 2). The acceleration of the bubble shell is reduced because of the radiative loss. The 
rise in the preshock sound speed and the acceleration of gas ahead of the shock reduced the 
shock compression (eqn. [4]). The shell becomes broader and less dense compared to the 
HD case (Fig. 4(c)). The sharp density jump at the shock front in the RHD case is also 
smoothed out before the gradient of e r diminishes to zero (Fig. 4(d)). The gas tends toward 
a freely expanding state in which the density profile is frozen into the flow. 

The diffusion effect speeds up the freeze-out of the ejecta structure. For our parameters, 
the freezing occurs at a velocity Vf ~ 1500 km s _1 and a time tf ~ 10 7 s, significantly earlier 
than the HD case where tf ~ 10 9 s and Vf ~ 1900 km s^ 1 (Fig. 4(f)). In the HD case, the 
expansion sustains until pressure equilibrium is reached. Both cases show a density contrast 
X > 10 between the shell and the preshock ejecta. The inner edge (contact discontinuity) of 
the shell is the densest with x ?t 100 . The shell thickness is increased to ~ 4% of the radius 
of the bubble. The swept-up mass in the shell remains about M s ~ 1 M (Fig. 4(g)). Note 
that we define the thickness as the distance from the contact discontinuity to the postshock 
outer dense region at the shock front, considering that the outermost postshock density 
profile flattens out into the ejecta substrate. At < 10 7 sec, the surface density of the shell 
drops below ~ 17 — 33g/cm 2 , the mean free path for the ~ lMeV 7-ray photons (Woosley, 
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Pinto, & Hartmann 1989; Shigeyama & Nomoto 1990). We stopped the RHD calculations 
prior to this time (Fig. 4(h)). When the initial Nickel mass is increased to 3 times 0.075 M & , 
the onset of freezing delays to tf ~ 10 8 s with v / ~ 1800 km s _1 and M s < 3 M & . For a 
10 times increase, tf ~ 10 8 s, Vf ~ 2400 km s _1 , and M s ~ 5 M Q (Fig. 5). It is therefore 
expected that the radiation relaxes before 10 8 s, given any large Ni overabundance. The 
reason is that the photon mean free path is increased by the broadening of the shell, and so 
diffusion becomes more rapid with time. 

3.2. Effect of Radiation 

The effect of radiative diffusion does not increase the density contrast of the Nickel 
bubble shell, as in shell formation during shock break-out in a Type II supernova where 
radiative energy losses by diffusion lead to a large compression that is eventually limited 
by the gas pressure (WC01, Chevalier 1976). In the latter case, the flow behind the shock 
wave is being decelerated, and the immediate postshock region is being compressed. A drop 
in the dominant radiation pressure then is compensated by compression of the gas to the 
point where gas pressure becomes important, resulting in a large density enhancement. In 
our case, the pressure gradient has an opposite sign; the pressure decreases outward from 
the contact discontinuity to the shock front, and the radiative loss to the preshock gas is 
replenished by the energy diffusing from the bubble. The diffusion of energy ahead of the 
shock front leads to a weaker shock front and less compression in the shell. 

3.3. Models for Type la supernovae 

We investigated the ejecta structure for Type la SNe using the inner flat component of 
the power law model, n = 8, E = 10 51 ergs, and M = 1.4 M . We present the evolution of 
the dynamical properties in Fig. 7. The simulations were our highest resolved runs initiated 
with t = 10 2 , Mni = 0.5 Mq, and uj = 1.5. The initiated time gives an initial mass density 
Po = 6.92 x 10 _1 g cm~ 3 , and an initial material and radiation energy density eo = e r o = 
5.5 x 10 14 erg cur 3 for a standard adiabat k — 6.1 x 10 14 (cgs units). Here a lower background 
pressure was used, e = 1.3 x 10 13 erg cm -3 . The initial 56 Ni has an accumulated radioactive 
energy density e ra do = 3.15 x 10 9 erg cm -3 . 

The evolution of the flow is similar to the case of core collapse supernovae. The Nickel 
bubble shell has a density increase toward the contact discontinuity, with the density contrast 
reaching x ~ 100. The diffusion effect speeds up the freeze-out of the ejecta structure. A 
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notable difference in the shell is that it is thinner, with a very prominent precurso (Fig. 6). 
Compared to core collapse supernovae, the ejecta has a higher ratio of background energy 
to mass (oc E/M), thus the acceleration of the bubble is weaker. The expansion rate a 
rises to a maximum 1.04 before radiation relaxes, barely faster than the free expansion. The 
thickness of the shell reaches (3 ~ 0.004 at the turnover. The freezing occurs at a velocity 
Vf ~ 7000 km s _1 and a time tf ~ 10 7 s. The gas becomes optically thin at ~ 10 6 sec. By 
this time, the RHD results do not diverge substantially from the HD results, but the velocity 
and the accumulated mass curves had started to flatten out, indicating that the radiation 
effect becomes dynamically less important. 

We note that the swept-up mass in the RHD case may be underestimated due to the 
large precurso; we expect M s ~ 0.5 M as in the HD case. For a less amount of 56 Ni, we 
found that M s < 0.3 M and v f ~ 5500 km s" 1 for M Ni = 0.225 M , and M s < 0.2 M 
and v f < 4000 km s" 1 for M Ni = 0.075 M (Fig. 5). 

3.4. Variation of the Initial Parameters 

We carried out simulations with a variety of initial conditions. We found the solu- 
tion in the saturated regime is insensitive to the initial background pressure in the ejecta 
(Fig. 8). In the case of a higher initial background pressure that exceeds the accumulated 
radioactive pressure, the sharp structure in the shell takes more time to build up (for which 
grid resolution also has an effect), but the solution gradually converges once the radioac- 
tive energy input becomes dominant. E.g., the case using a smaller background pressure 
P = 2.7 x 10 10 dyn cm -2 is no different from the standard case using P = 2.7 x 10 18 dyn cm~ 2 
(for the core collapse SN model with t — 10 2 sec), except that the former shows a more 
pronounced postshock oscillation trailing in the bubble. Likewise, the solution in our Type 
la model appears insensitive to the variation of background pressure. Our case using the 
standard k shows little differences from the case in Fig. 7. We expect that the insensitivity 
would sustain with a higher value of k. Simulations initialized before 10 4 sec allow the flow 
to saturate before the age of interest, and so the results are also insensitive to the initial dis- 
tribution of energy between radiation and material, absorption rate, and details of radiative 
transport. 

We also varied the initial 56 Ni density contrast. We assumed 56 Ni initially in pressure 
equilibrium with the ambient ejecta, so that its initial background energy increases with 
decreasing u. Thus for ou — 0.1 the shell expands to an outer radius. However, the change 
from uj = 3 to u = 1 does not alter the shell properties we determined (Fig. 8). 
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4. ANALOGY TO PULSAR WIND NEBULAE 

Chevalier (1984) has considered a constant-luminosity pulsar bubble interacting with 
constant density, freely-expanding supernova ejecta. We compare his self-similar solution 
with our hydrodynamical case, since the time evolution of the driving power is the only 
difference. The self-similarity exists for a power-law ejecta density p oc v~ n t~ 3 and a time- 
varying pulsar luminosity L oc t~ l , when n < 3 and I < 1 (Jun 1998); the expansion law is 
r oc t( 6 - n -0/(5-n) ) ag obtain^ from dimensional analysis. For a constant luminosity I = 
and flat ejecta n = 0, r oc t 12 . The shock cannot be slowed down in the ejecta more than 
the free expansion rate, so the self-similar solutions only exist for a > 1, i.e. I < 1. 

In our solutions, the shell starts to accelerate with a > 1. The evolution of the density 
profile can thus be compared to the self-similar solutions with varying /. In the self-similar 
solutions, decreasing I (i.e. increasing power input), or increasing the expansion rate, in- 
creases the shell thickness. For / = and a = 1.2, we have (3 = 0.02. The shell thickness in 
our simulations shows a more complicated evolution than indicated by the measured expan- 
sion rate; the maximum expansion rate measured occurs at 6 x 10 5 sec, while the maximum 
thickness occurs at 5 x 10 6 sec. The first turnover in the thickness evolution seems to cor- 
relate with the change in the luminosity power index (see Fig. 9) at 5 x 10 6 sec; the second 
turnover occurs after 10 7 sec when the self-similarity breaks down because I > 1. 

We overplot our solution at 4 x 10 5 sec on the self-similar solution with I = 0.03 in 
Fig. 10. The density distribution in the shell shows a sharp inward increase towards the 
contact discontinuity, with over 80% of the shell mass concentrated within the inner 50% 
region in radius (Fig. 11). The self-similar density is infinite at the contact discontinuity. 
The coarseness of the grid obviously sets limits on the highest density computed in the grid 
domain. Our results show a higher density contrast on fixed grids, because the numerical 
noise is smaller than on expanding grids. The shell can be distorted by instabilities (Basko 
1994; Jun 1998), but we expect the density contrast to remain comparable. 

Basko's (1994, his Fig. 1) HD solutions suggest that the density contrast of the shell 
is only x ~ 5 at 4 x 10 5 sec. We reproduced his case and found that the shell broadens to 
(3 > 0.1 at ~ 4 x 10 4 s, the turnover is much earlier than in our calculation using to — 10 2 
s. In his case using t = 10 4 s, the radioactive energy accumulated to 10 4 s was deposited 
all at one time in the beginning of the calculation, so an excessive acceleration and a large 
broadening of the shell were the result. Basko's Fig. 6 for the shell surface density has the 
same earlier turnover, suggesting this heating effect. 

We note that the shell characteristics such as the density contrast and the thickness 
ratio should only be weakly dependent on the initial 56 Ni mass, since these properties can 
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be approximated by the self-similar solutions. 

The deformation of the bubble-shell interface due to the Rayleigh- Taylor instability is 
likely to have a morphology similar to that of the pulsar bubble shell studied by Jun (1998) 
in two-dimensional simulations. The main difference is that magnetic fields constrain the 
instability in the pulsar bubble case, but are unlikely to be important for our case. The 
Rayleigh- Taylor spikes are directed inward, because of the shell acceleration. However, the 
instability should eventually end when the shell stops accelerating in the comoving stage. 
The effect of radiation diffusion is to smooth the pressure gradients, which ends the shell 
acceleration and the instability growth. The Jun (1998) simulations indicate that the insta- 
bility does not cause radioactive material to have a significantly higher final velocity (see 
also Basko 1994). 



5. NICKEL BUBBLE SHELL AS EJECTA CLUMPS 

We now compare the properties of the Ni bubble shell with the initial properties of the 
clumps/bullets inferred from our previous clump-remnant simulations for Tycho's remnant 
and the Vela supernova remnant. We consider that the ejecta clumps originate as components 
from the breakup of the shell, and examine if the inferred clumps (shell fragments) are able 
to survive the crushing in the remnant. The robustness of the clumps can be determined by 
three parameters: the initial density contrast, x, between the clump and supernova ejecta; 
the initial impact time with the reverse shock, related to the initial ejection velocity; and 
the initial size of the clump. The high density contrast created across the shell, x ~ 100) 
is compatible with that of ejecta clumps needed to survive crushing in the remnant. In our 
model for the Vela remnant (M = 10 M , E 51 = 1, and M Ni = 0.075 M ), the shell's 
frozen-in velocity, Vf ~ 1500 km s -1 , is within the range 1000 < v < 3000 — 4000 km s" 1 we 
determined as the most likely origin for the Vela bullets (WC02). If the shell has not been 
greatly disturbed, the clumps should be present in bands at the restricted velocity space; 
they are expected to enter the remnant's intershock region at a normalized time t' ~ 2.2, 
or ~ 2800 yrs (see Fig. 1 in WC02 and Fig. 12). The velocity-restricted nature may have 
been detected in the wavelike clumpy structure immediately below Tycho's blast wave front 
(Reynoso et al. 1999, WC01). If the clumps attain a size about the thickness of the shell, 
(3 = 0.04, then at the time of the initial reverse shock impact, the clumps have an initial size 
do ~ 8% m terms of the intershock width, or r < 3% in terms of the forward shock radius 
(Fig. 1 in WC02). Compared to the late clump- remnant interaction in WC02 that used a 
larger clump size do ~ 1/3, the inferred clumps do not appear to have a sufficient strength to 
cause a protrusion on the forward shock. Nevertheless, since the computed highest density 
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is limited by numerical resolution, a larger density contrast should be achieved in a small 
region. An increase in the density contrast, e.g., from 100 to 1000, can significantly improve 

— 1/2 

the strength, considering that the cloud crushing time scale is proportional to t cc oc r c 
for a fixed mass (WC02, Klein et al. 1994). Such small-sized but dense clumps are likely to 
form, given the steep density profile of the shell. 

In our model for Type la supernovae (M = 1.4 M , E 51 = 1, and M Ni = 0.5 M Q ), 
the frozen-in velocity reaches Vf ~ 7000 km s _1 , resulting in a clump-remnant interaction 
initiated in an early state, t' ~ 1.02. Compared to the present state t' ~ 1.7 of Tycho's 
remnant, the ejection has a strength as indicated for Tycho's knots in our previous 2-D HD 
simulations (WC01). In addition, the swept-up mass ~ 0.5 M in the shell is reasonably 
in excess of the mass 0.002 M for the Si+S clump and 0.0004 M for the Fe clump in 
Tycho's SE outline (Hwang, Hughes, & Petre 1998). However, there is evidence for Tycho 
that it is subluminous (van den Bergh 1993; Schaefer 1996), implying a deficiency in 56 Ni. 
Observations of SN 1006 suggest an initial Nickel mass in the range of 0.075 — 0.16 M 
(Hamilton et al. 1997, Mazzali et al. 1997). For Sn 1991 bg, Mazzali et al. (1997) estimate 
0.07 Mq of 56 Ni extending to a velocity of 5,000 km s -1 . In our model, the velocity only 
extends to < 4, 000 km s" 1 with M Ni = 0.07 M . 

For core collapse supernovae, we consider that the object RX J0852-4622 superposed the 
boundary of the Vela remnant is a result of the clump-remnant interaction. WC02 estimated 
a mass of 0.1/3 M , where j3 < 1, for the responsible clump, which gives a fraction < 10% in 
terms of the total swept-up mass in the shell. This mass fraction appears somewhat large 
compared to the Si+S clump in Tycho's remnant. It is unclear how such a massive clump is 
created, if the initial 56 Ni is not overabundant. 

In our scenario, the initial 56 Ni has a symmetric distribution in the center of supernova 
ejecta. It should be noted that dynamical processes prior to the Nickel bubble effect, e.g., 
neutrino-driven convection (Kifonidis et al. 2000), may give rise to an anisotropic distribu- 
tion of elements. If 56 Ni is mixed outwards during the explosion, the amount of radioactive 
material inside each individual bubble could be orders of magnitude lower; the subsequent 
radioactive heating may cause merging between Nickel bubbles, and so give further com- 
pression. This process remains to be examined by future multidimensional calculations. 
Our result here should be noted as a lower limit on the ejecta-clump properties for the 
clump-remnant interaction. 
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6. CONCLUSIONS 

We investigated the ejecta structure in supernovae due to the radioactive heating from 
the 56 Ni — > 56 Co — > 56 Fe decay sequence, under the assumption of spherically symmetric 
flow. Two approaches were used: hydrodyamics (HD) simulations and radiation hydrody- 
namics (RHD) simulations. The RHD calculation included a proper treatment of radiation 
transport prior to the age 10 7 sec, when about 25% of the total radioactive energy is still 
to be released. Radiative effect becomes dynamically less important as the evolution tends 
towards the optically thin regime. Based on the mutual coherence in the HD and RHD 
calculations, we obtained the properties of the shell driven by the radioative expansion. We 
believe that going to later stages of diffusion would not change this conclusion. 

The expansion of the Nickel bubble can sweep up a dense shell of < 1 M shocked 
ambient ejecta in core collapse and Type la supernovae, with the highest density in the shell 
over 100 times the ambient ejecta density. The shell has an inward density increase toward 
the bubble-shell interface, and the highest density contrast computed is limited by numerical 
resolution. For the hydrodynamical solution, the structure of the shell can be approximately 
described by a self-similar solution. For the RHD solution, the main difference is that 
radiative diffusion gives rise to a broader and less dense shell that freezes out in the ejecta. 
Because the radioactive pressure eventually dominates, the saturated solution is not sensitive 
to the initial conditions such as the background pressure and the initial density contrast of 
Ni; the only crucial factor to the solution is the initial 56 Ni abundance. We presume that 
a higher background pressure due to the uncertainty of the adiabat would not affect the 
results. 

We examined if the properties of the shell fragments are comparable to that of the 
ejecta clumps protruding the outlines of Tycho's remnant and the Vela remnant. The density 
contrast created across the shell, x ~ 100, is compatible with that of the ejecta clumps needed 
to survive the clump-remnant interaction. In our standard Type la model, the shell's frozen- 
in velocity attained at ~ 10 6 sec appears sufficient to give an ejection as indicated for Tycho's 
knots, in spite that most (< 75%) of the radioactive energy is still to be released. For core 
collapse supernovae like the Vela remnant, however, the small size and the late ejection of 
the inferred clump result in a weak strength to resist crushing and further cause protrusions 
on the remnant outline. 

I am grateful to Roger Chevalier for supporting this work, and Toshi Shigeyama for 
useful correspondence on Type la supernovae. I thank the referee for careful and constructive 
comments on the manuscript. The computations were carried out on the IBM SP2 at 
University of Virginia, and the PC cluster provided by Taiwan's NSC grant NSC93-2112- 
M033-002. Support for this work was also provided by NSC grant NSC93-2112-M033-013.. 
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Fig. 1. — Evolution of the instantaneous and accumulated heating rate of the Nickel and Cobalt 
radioactive decay. 
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Fig. 2. — Distributions of density, velocity, gas energy density, and radiation energy density of the 
Nickel bubble structure in the flat-density ejecta of the core collapse supernova at 10 5 and 10 6 sec. 
The initial Nickel mass is M n = 0.075 M & . The grid has 8000 uniform zones that resolve the shell 
into ~ 100 zones. 
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Fig. 3. — Density profiles of the Nickel bubble shell in the core collapse SN model at 10 5 and 10 6 
sec. a: HD solution on a nonuniform grid of 8000 zones that resolves the shell into 200 zones, b: 
HD solution on a uniform grid of 8000 zones, c: RHD solution on a uniform grid of 8000 zones, d: 
RHD solution on a uniform grid of 2000 zones. 
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Fig. 4. — Evolution of the dynamical properties of the Nickel bubble shell in the core collapse SN 
model. Solid and dotted lines represent the RHD and the HD case, respectively, (a) Radius of the 
shock R s h- (b) Expansion rate of the shell, a = dlnR s h/dlnt, equivalent to the velocity contrast 
between the shock and the ambient free-expanding ejecta. (c) Thickness ratio (3 of the shell. The 
value is influenced by numerical resolution, (d) Density contrast between the densest region in the 
shell and the ejecta substrate. Dashed and dash dot lines represent the average density contrast in 
the HD and the RHD case, respectively, (e) Density contrast between the bubble and the ejecta 
subtrate. (f) Shell velocity v s h = dR s h/dt (left) and flow velocity Vfi ow = R s h/t (right) at the shock 
front, (g) Swept-up mass M s in the shell. In the RHD case, two estimates are given for M s due to 
the lower resolution. Solid line: estimated within the inner sharp edge of the shock front. Dashed 
line: estimated within the outer edge of the shock front where the density flattens out. (h) Surface 
density a s = M s jAmR 2 sh of the shell. 
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Fig. 5. — Evolution of shell and shock velocity and accumulated mass in the shell with varying 
initial Ni abundance. Solid and dotted lines represent the RHD and the HD case, respectively. 
Left: = 0.075 (bottom line), 0.225 (middle line), and 0.75 M Q (top line) for the core collapse 
SN model. Right: M Ni = 0.075 (bottom line), 0.225, (middle line), and 0.5 M Q (top line) for the 
Type la SN model. 
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Fig. 6. — Distributions of density, velocity, gas energy density, and radiation energy density of the 
Nickel bubble structure in the Type la SN model at 10 5 and 10 6 sec. The initial Nickel mass is 
M n = 0.5 Mq. The grid has 16000 uniform zones. 
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Fig. 7. — Evolution of the properties of the Nickel bubble shell in the Type la SN model. Solid 
and dotted lines represent the RHD and the HD case, respectively, (a) Radius of the shock R s h- 
(b) Expansion rate of the shell, a = dlnR s h/dlnt, equivalent to the velocity contrast between the 
shock and the ambient free-expanding ejecta. (c) Thickness ratio (3 of the shell. The thickness is 
estimated within the densest point and the outer edge of the shock front. In the RHD case it is 
not well defined due to a low resolution, (d) Density contrast between the densest region in the 
shell and the ejecta substrate. Dashed and dash dot lines represent the average density contrast in 
the HD and the RHD case, respectively, (e) Density contrast between the bubble and the ejecta 
subtrate. (f) Shell velocity v s h = dR s h/dt (left) and flow velocity Vfi ow = R s h/t (right) at the shock 
front, (g) Swept-up mass M s in the shell. In the RHD case, M s is estimated within the outer edge 
of the shock front without the precurso. (h) Surface density a s = M s /ATrR 2 sh of the shell. 
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Fig. 8. — Density profiles at 10 6 sec (top) and evolution of shell and flow velocities (bottom) in 
HD models for core collapse SNe with varying initial parameters. The runs in dotted and dash dot 
lines use a resolution of 8000 zones, others use a resolution of 2000 zones. Solid line: to = 10 4 sec, 
e = 2.7 x 10 11 (cgs units), and uj = 3. Dotted line: t = 10 2 
uj = 3. Dashed line: to = 10 4 sec, eo = 2.7 x 10 11 (cgs units), 
eo = 2.7 x 10 18 (cgs units), and uj = 1. Dash dot dot dot: to = 
and uj = 0.1. 



sec, eo = 2.7 x 10 18 (cgs units), and 
and uj = 1. Dash dot: to = 10 2 sec, 
-- 10 2 sec, e = 2.7 x 10 18 (cgs units), 
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Fig. 9. — Luminosity index I of Nickel radioactivity and expansion rate a derived from the self- 
similar solution, a = (6 — n — /)/(5 — n). 




Fig. 10. — Hydrodynamical solutions of the Nickel bubble shell in the core collapse SN model at 
4.0 x 10 5 sec overplotted on the self-similar solutions of 7 = 4/3, n = 0, I = 0.3, and a = 1.14. The 
velocity and pressure show large post shock oscillations behind the contact discontinuity. 
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Fig. 11. — Top: Density contrast of the self-similar solution at 4 x 10 5 sec integrated over a 
radial fraction of the shell from the contact discontinuity. The innermost 1% region of the shell 
has an average density contrast over 600, while the total average is 20. Bottom: Integrated mass 
distribution. 
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Fig. 12. — Evolution of velocities in a supernova remnant with n=8 ejecta and a constant density- 
ambient medium as described in WC02. Dotted line: flow velocity of the ejecta immediately below 
the reverse shock (RS). Solid lines: the pattern velocities and deceleration parameters 5 = dlnr/dlnt 
(inset) of the forward shock (FS), the contact discontinuity (CD), and the reverse shock with time. 
The quantities v' and t' are normalized to scaling parameters given in WC 02. In our model 
for core collapse supernovae, v = 3162 v' (km/s) and t = 1271 t' (yr); for Type la supernovae, 
v = 8452 v' (km/s) and t = 244 t' (yr). 
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